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Abstract. Transition probabilities governing the inter- 
action of energy packets and matter are derived that al- 
low Monte Carlo NLTE transfer codes to be constructed 
without simplifying the treatment of line formation. These 
probabilities are such that the Monte Carlo calculation 
asymptotically recovers the local emissivity of a gas in 
statistical equilibrium. Numerical experiments with one- 
point statistical equilibrium problems for Fe II and Hy- 
drogen confirm this asymptotic behaviour. In addition, 
the resulting Monte Carlo emissivities are shown to be far 
less sensitive to errors in the populations of the emitting 
levels than are the values obtained with the basic emissiv- 
ity formula. 

Key words: methods: numerical - radiative transfer - 
stars: atmospheres - Line: formation 



<* 1. Introduction 



* 



When Monte Carlo methods are used to compute the spec- 
tra of astronomical sources, it is advantageous to work 
with indivisible monochromatic packets of radiant energy 
and to impose the constraint that, when interacting with 
matter, their energy is conserved in the co-moving frame. 
The first of these constraints leads to simple code and the 
second facilitates convergence to an accurate temperature 
stratification. 

For a static atmosphere, the energy-conservation con- 
straint automatically gives a divergence- free radiative flux 
even when the temperature stratification differs from the 
radiative equilibrium solution. A remarkable consequence 
is that the simple A-iteration device of adjusting the tem- 
perature to bring the matter into thermal equilibrium with 
the Monte Carlo radiation field results in rapid conver- 
gence to the close neighbourhood of the radiative equilib- 
rium solution (Lucy 1999a). An especially notable aspect 
of this success is that this temperature-correction proce- 
dure is geometry-independent, and so these methods read- 
ily generalize to 2- and 3-D problems. 
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For an atmosphere in differential motion, the energy- 
conservation constraint yields a radiative flux that is rig- 
orously divergence- free in every local matter frame. Deter- 
mining the temperature stratification by bringing matter 
into thermal equilibrium with such a radiation field - i.e., 
by imposing radiative equilibrium in the co-moving frame 
- is an excellent approximation if the local cooling time 
scale is short compared to the local expansion time scale. 
This condition is well satisfied for the spectrum-forming 
layers of supernovae (SNe) and of hot star winds (Klein & 
Castor 1978). 

The constraint that the energy packets be indivisible is 
advantageous from the point of view of coding simplicity. 
The interaction histories of the packets are then followed 
one-by-one as they propagate through the computational 
domain, with there being no necessity to return to any of a 
packet's interactions in order to continue or complete that 
interaction. This is to be contrasted with a Monte Carlo 
code that directly simulates physical processes by taking 
its quanta to be a sampling of the individual photons. Ab- 
sorption of a Monte Carlo quantum is then often followed 
by the emission of several quanta as an atom cascades 
back to its ground state. Multiple returns to this interac- 
tion are then necessary in order to follow the subsequent 
paths of each of these cascade quanta. The resulting cod- 
ing complexity is of course compounded by some of these 
quanta creating further cascades. 

Although coding simplicity argues strongly for indi- 
visible packets, a counter argument is the apparent im- 
plied need to approximate the treatment of line formation. 
Thus, in Monte Carlo codes for studying the dynamics of 
stellar winds (Abbott & Lucy 1985; Lucy & Abbott 1993) 
or for synthesizing the spectra of SNe (Lucy 1987; Mazzali 
& Lucy 1993), the integrity of the packets could readily 
be maintained since lines were assumed to form by co- 
herent scattering in the matter frame. But significantly, 
an improved SN code has recently been described (Lucy 
1999b) in which branching into the alternative downward 
transitions is properly taken into account without sacri- 
ficing indivisibility. Accordingly, an obvious question now 
is whether Monte Carlo techniques can be developed that 
enforce energy-packet indivisibility and yet do not have to 
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adopt any simplifications with regard to line formation. If 
this can be achieved, then Monte Carlo codes for general 
NLTE transfer problems become feasible. 

2. Macro-atoms 

As discussed in Sect. 1, it is common in Monte Carlo 
transfer codes to quantize radiation into monochromatic 
energy packets. But matter is not quantized, neither nat- 
urally into individual atoms nor artificially into parcels of 
matter. Instead, the continuum description of matter is 
retained, with macroscopic absorption and scattering co- 
efficients governing the interaction histories of the energy 
packets. 

Nevertheless, it now proves useful to imagine that mat- 
ter is quantized into macro-atoms whose properties are 
such that their interactions with energy packets asymp- 
totically reproduce the emissivity of a gas in statistical 
equilibrium. But these macro-atoms, unlike energy pack- 
ets, do not explicitly appear in the Monte Carlo code. As 
conceptual constructs, they facilitate the derivation and 
implementation of the Monte Carlo transition probabili- 
ties that allow in an accurate treatment of line formation. 

The general properties of macro-atoms are as follows: 

1) Each macro-atom has discrete internal states in one- 
to-one correspondence with the energy levels of the atomic 
species being represented. 

2) An inactive macro-atom can be activated to one of 
its internal states i by absorbing a packet of kinetic energy 
or a packet of radiant energy of an appropriate co-moving 
frequency. 

3) An active macro-atom can undergo an internal tran- 
sition from state i to any other state j without absorbing 
or emitting an energy packet. 

4) An active macro-atom becomes inactive by emitting 
a packet of kinetic energy or a packet of radiant energy of 
an appropriate co-moving frequency. 

5) The de-activating packet has the same energy in the 
macro- atom's frame as the original activating packet. 

Figure 1 illustrates these general rules. An inactive 
macro-atom, with internal states shown schematically, en- 
counters a packet of energy eo and is activated to one of 
these states. The active macro-atom then undergoes two 
internal transitions before de-activating itself by emitting 
a packet of energy e - 

Subsequently, energy packets will in general be referred 
to as e-packets but also as r- or fc-packets when specifying 
their contents to be radiant or kinetic energy, respectively. 

3. Transition probabilities 

In Sect. 2, the concept of a macro-atom was introduced 
by stating some general properties concerning its inter- 
action with e-packcts. The challenge now is to derive ex- 
plicit rules governing a macro-atom's activation, its subse- 
quent internal transitions, and its eventual de-activation. 
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Fig. 1. Schematic representation of the interaction of a 
macro-atom with a packet of energy eo- The macro atom 
is activated by absorbing the energy packet, makes two 
internal transitions, and then de-activates by emitting a 
packet of energy eo . 



Asymptotically, the result of obeying these rules must be 
the emissivity corresponding to statistical equilibrium. 



3.1. Energy flow rates 

For the moment, we drop the notion of a macro-atom and 
consider a real atomic species interacting with its environ- 
ment. Let £j denote the excitation plus ionization energy of 
level i and let Rij denote the radiative rate for the transi- 
tion i — » j. The rates per unit volume at which transitions 
into and out of i absorb and emit radiant energy are then 



A? 



Rata and E, — Rueu 



(1) 



respectively, where ea — hvn = Ci — eg. Note the sum- 
mation convention adopted for the suffix £, which ranges 
over all levels < i, including those of lower ions. Similarly, 
below, the suffix u implies summation over all levels > i, 
including those of higher ions. 

The corresponding rates at which kinetic energy is ab- 
sorbed from, or contributed to, the thermal pool by tran- 
sitions to and from level i are 



A? 



Caei 



and 



E? 



Cue 



il^il j 



(2) 



where CV, is the collisional rate per unit volume for the 
transition i — > j. 

If we now define the total rate for the transition i — ► j 
to be IZij = Rij + Cij , then the net rate at which level i 
absorbs energy is 



if + A? - E* - E? = (1Z U - TluWi - u) 



(3) 



This is an identity that follows directly from the defining 
Eqs.(l) and (2); it is therefore quite general and does not 
assume statistical equilibrium. 
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3.2. Statistical equilibrium 

We now assume that the level populations m are in sta- 
tistical equilibrium. For level i, this implies that 



(Ku -Hu) + {K m -K lu ) = Q 



(4) 



A useful alternative representation of statistical equi- 
librium is obtained by multiplying Eq.(4) by e, and then 
eliminating the term (TZu—TZu)^ using Eq.(3). The result 
can be written in the form 



rlt 



?C 



Af + Af + TZuiEi + TZ u e e . 



(5) 



Eq.(4), the conventional equation of statistical equilib- 
rium, balances the rates at which basic atomic processes 
excite and de-excite level i. As such, it directly relates 
to Nature's quantization of radiation into photons and of 
matter into atoms. In contrast, Eq.(5), though mathemat- 
ically equivalent, deals with macroscopic energy flow rates 
in a finite volume element. These flows can now be quan- 
tized into indivisible e-packets. Moreover, we can think of 
the volume element as a macro-atom with discrete energy 
states. 



3.3. Interpretation 

Eq.(5) expresses the fact that in statistical equilibrium 
the contribution from level % to the energy content of unit 
volume is stationary. In consequence, the net rate at which 
level % gains energy - the right-hand side of Eq.(5) - equals 
the net rate of loss - the left-hand side. 

But the importance here of Eq.(5) lies in the various 
terms contributing to gains and losses by level i and their 
relevance for constucting transition rules for macro-atoms. 
The net rate of gain comprises the expected absorption 
terms Af- and Af plus the terms !Z U i£i and TZaei that 
clearly represent energy flows into i from upper and lower 
levels. Similarly, the net rate of loss comprises the ex- 
pected emission terms E^ and Ef plus the terms TZi U €i 
and IZuei representing energy flows out of % to upper and 
lower levels. 

The above remarks imply definitive values for the en- 
ergy flows between level i and other levels. But this is not 
true. If Eq.(4) is rewritten as 



Ti-iu + r^il — TZ-ui + 'R'l 



(6) 



then comparison with Eq.(5) shows immediately that an 
arbitrary quantity of energy e may be added to ej and ££ 
without invalidating this equation. But this merely shifts 
the zero point of the energy scale for excitation and ion- 
ization, which we are always free to do. Nevertheless, this 
freedom implies a corresponding indefinitcness in the en- 
ergy flow rates between levels. 



3.4- Stochastic interpretation 

Notwithstanding this indefinitcness, we now interpret 
Eq. (5) in terms of macro-atoms absorbing and emitting e- 
packets or undergoing transitions between internal states. 
In this interpretation, the terms Af- and Af obviously 
represent the activation of macro-atoms to state i due to 
the absorption of r-packets and of fc-packets, respectively. 
Now consider an ensemble of active macro-atoms in 
state i. For this ensemble to reproduce the behaviour of 
the real system, the relative numbers of the macro-atoms 
that subsequently de-activate with the emission an r- or 
fc-packet or which make a transition to another internal 
state must be proportional to the relative values of the 
corresponding terms on the left-hand side of Eq.(5). Ac- 
cordingly, for an individual macro-atom in state i, the 
probabilities that it de-activates with the emission of an 
r-packet or a fc-packet are 



Pi 



" - E?/Di and pf = E? / D t , 



where 
Di = E? 



E 



he 



- TZ-iuCi + Tlii^f. = {Ti-u + TZi U )ei 



(7) 



(8) 



Similarly, the probabilities that it makes an internal tran- 
sition to particular upper or lower states are 



Pi 



TZiuSi/Di and p ie = 7£^q/A 



(9) 



Unlike transition probabilities for real atoms, these 
analogues for macro- atoms depend on ambient conditions. 
Consequently, in the course of a NLTE calculation, they 
are iterated on just as are Eddington factors in various 
other radiative tranfer schemes (Auer & Mihalas 1970; 
Hummer & Rybicki 1971). Moreover, as with Eddington 
factors, the Monte Carlo transition probabilities are di- 
mcnsionless ratios that are likely to converge faster than 
do their dimensional numerators and denominators. 

3.5. Excitation equilibrium 

When Eq.(5) is summed over all energy levels, the energy 
flows between different levels cancel, giving 



J2(A« + Af) = J2(E« + E?) 



(10) 



Thus, in statistical equilibrium, the energy stored in the 
form of excitation and ionization is stationary. For the 
macro-atoms, this is obeyed rigorously by each activation 
- de-activation event since the emitted packet's energy 
equals that of the absorbed packet - see Figure 1. 

4. Alternative formulations 

Monte Carlo transition probabilities have been defined in 
Sect. 3, but their non-negativity was not established. Of 
concern in this regard is stimulated emission when level 
populations are inverted. However, in anticipation of this 
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issue, radiative rates were introduced without specifying 
whether stimulated emission contributes positively to Rij 
or negatively to Rji. We now exploit this flexibility in 
order to avoid negative probabilities. 



4-1- General case 

In the general case, inverted level populations may occur 
- i.e., gjUi > gifij for some i > j . 

4.1.1. Definitions of rates 

In order to prevent the probabilities becoming negative 
when levels invert, stimulated emissions must be added to 
spontaneous emissions and not treated as negative absorp- 
tions. Accordingly, for bound-bound (b-b) transitions, the 
radiative rates per unit volume arc defined to be 



Rij = (Aij + BijJfAm and Rji = B^J^nj 



(11) 



where J? and </"■ are the mean intensities averaged over 
the line's emission and absorption profiles - see Mihalas 
(1978, p78). Similarly, for free-bound (f-b) and bound-free 
(b-f) transitions, we define 



R* 



- („, S P 



(a t p + af )n K n e and R iK = Jim 



(12) 



Here af and af are the rate coefficients for spontaneous 
and stimulated recombinations to level i, and 7$ is the un- 
corrected rate coefficient for photoionizations from level 
i. Each of these three quantities can be expressed as an 
integral over frequency involving the b-f absorption coef- 
ficient for an atom excited to level i - see Mihalas (1978, 
ppl30-131). 

For collisions, a population inversion gives a negative 
rate if de-excitations are treated as negative excitations. 
This is avoided by defining 



y-^ij — Qijnin e 



and 



^ji — Qjinjn e 



(13) 



With these expressions for the radiative and collisional 
rates, the probabilities defined by Eqs.(7) and (9) are non- 
negative provided only that the e^'s are non-negative. This 
latter condition is satisfied with the standard convention 
that the ground state of the neutral atom has zero exci- 
tation energy. 



4.1.3. Emission of packets 

If the Monte Carlo transition probabilities result in a 
macro-atom de-activating radiatively from state i, the 
next step is to determine the frequency of the photons 
comprising the emitted r-packet. First we suppose that % 
corresponds to a bound level. 

Because Ra and therefore Ef here include stimu- 
lated emission, the process that radiatively de-activates 
the macro-atom may be either a spontaneous or a stim- 
ulated emission. The ratio of the probabilities of these 
alternatives is q = Ef /Ef, where 



El 



Aunitu 



Ef and Ef 



BuJu n i e U 



771st 



(14) 



are the contributions to Ef- from spontaneous and stim- 
ulated emissions. Knowing q, we can choose between the 
two alternatives with a standard Monte Carlo procedure. 
Thus, if x is a random number from the interval (0, 1), we 
select spontaneous emission if a; < q/(l + q) and stimu- 
lated otherwise. 

Having thus decided the emission process, we must 
next choose a downward transition. For spontaneous line 
emission, the transition i — > j is selected with probability 
Ef /Ef . For stimulated emission, on the other hand, the 
selection probability is Ef: /Ef. 

With the transition thus determined, the frequency v 
of the r-packet is selected by randomly sampling the line's 
emission profile </%. Thus, if x again denotes a random 
number from (0, 1), then v is determined by the equation 



(jf v dv = x 



(15) 



This equation can of course always be solved numerically 
for v. However, elegant and efficient procedures for sam- 
pling standard profiles are available (Lee 1974a, b). 

Now we consider a macro-atom that de-activates from 
a continuum state k. In this case, the probabilities of spon- 
taneous and stimulated emission are in the ratio Ef : Ef , 
where 



Ef 



sp 



H n K e K £ 



and 



& K = a u n K e K t 



(16) 



are the contributions to Ef from spontaneous and stim- 
ulated emissions. Thus v is selected by first deciding be- 
tween spontaneous and stimulated emission and then ran- 
domly sampling the energy distribution of the chosen pro- 
cess's recombination continua. 



4.1.2. Absorption of packets 

Because Ru and therefore Af are here defined without 
correcting for stimulated emission, the macroscopic line- 
and continuum-absorption coefficients that determine the 
flight paths of r-packets must also be defined without this 
correction. This ensures a positive absorption coefficient 
even for a transition with a population inversion. 



4.1.4. Direction of propagation 

If the above selection procedure rules that an r-packet is 
emitted spontaneously, then a new direction of propaga- 
tion is chosen in accordance with this process's isotropic 
emission. On the other hand, for stimulated emission, the 
new direction of propagation is that of the stimulating 
photon. Thus, the new direction will be in solid angle 
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duo at 0,cj> with probability duj/Air x l v (0,<j>)/ J v , where 
v is the frequency of the emitted r-packet. Accordingly, 
a Monte Carlo code that treats stimulated emission sepa- 
rately must store a complete description of the radiation 
field -i.e., I v {0,4>). 

4-2. Standard case 

For problems where population inversions are not antic- 
ipated, we can usefully make the traditional assumption 
that lines have identical emission and absorption profiles 
and treat stimulated emissions as negative absorptions - 
sec Mihalas (1978, p78). 

4.2.1. Definitions of rates 

The radiative rates for b-b transitions are then 

Rij = AijHi and R ]% = (BjiUj - B lj n i )J : j l . (17) 

Similarly, for f-b and b-f transitions, we define 



Rk 



sp 



a^n^e and R iK = 7? 



(18) 



where the photionization coefficient is now corrected for 
stimulated recombinations. 

For collisions, the absence of population inversions 
allows us to treat de-excitations as negative excitations 
without the risk that Eqs.(7) and (9) will give negative 
probabilities. Accordingly, we now define 



Cij = and Cji = {qjifij - qijHiji 



(19) 



This then implies that Ef and therefore also pf = for all 
i. Energy transfer from the radiation field to the thermal 
pool then occurs explicitly only via f-f absorptions. 

4.2.2. Absorption of packets 

Because Ri u and therefore Af are here defined with the 
correction for stimulated emission included, the macro- 
scopic line- and continuum-absorption coefficients must 
also include this correction. In the posited absence of pop- 
ulation inversions, these absorption coefficients are posi- 
tive. 

4.2.3. Emission of packets 

Because Ru and therefore Ef- now exclude stimulated 
emission, the process that radiatively de-activates a 
macro-atom is always a spontaneous emission. If i is a 
bound state, the frequency v of the emitted r-packet is 
then decided as follows: the transition i — > j is selected 
with probability A^niCij / 'E^ , and then v is selected by 
randomly sampling this transition's emission profile, as in 
Sect. 4.1.3. 

For de-activation from a continuum state, v is selected 
by randomly sampling the energy distribution of the spon- 
taneous recombination continua. 



4.2.4. Direction of propagation 

Because the de-activating process is in this case sponta- 
neous emission, the new direction of propagation is se- 
lected according to isotropic emission. Thus, we now do 
not need to store l u (0,(f>). In fact, from the Monte Carlo 
radiation field generated at one iteration, we only require 
the mean intensities J„ . These allow us to compute transi- 
tion probabilities from Eqs.(7) and (9) for use during the 
next iteration. 



4-3. Large velocity gradients 

The procedures described in Sects. 4.1 and 4.2 apply to 
both static and moving media. But for some important 
problems involving moving media, a substantial speeding 
up of the calculation with negligible loss of accuracy is 
possible by applying Sobolev's theory of line formation. 
In doing so, we take advantage of a small dimensionlcss 
quantity - the ratio of a line's Doppler width to the typical 
flow velocity, which implies an essentially constant velocity 
gradient over the zone in which a given pencil of radiation 
interacts with a particular line. The Monte Carlo codes for 
hot star winds and SNe cited in Sect. 1 treat line formation 
in the Sobolev approximation. 

The simplest case of this kind is that of homologous 
spherical expansion, as is commonly assumed for SNe. 
This case will be treated here since it will be used in 
the test calculations of Sect. 5. But generalization to a 
spherically-symmetric stellar wind is readily carried out 
by referring to Castor & Klein (1978). We also assume no 
population inversions and so treat stimulated emissions as 
negative absorptions, as in Sect. 4.2. 

4.3.1. Definitions of rates 

The radiative rates for b-b transitions arc then 



Rj 



AijPjifii and R jt = (B 3i n 3 - Bijn^PjiJ^ . (20) 



Here J^ is the mean intensity at the far blue wing of the 
transition j — > i, and (3ji is the Sobolev escape probability 
for this transition, given by 



Tji 

where Tji, the transition's Sobolev optical depth, is 



iji 



{BjiHj - B i3 n 



l 3 >H j 



hctE 

47T 



(21) 



(22) 



with t E being the elapsed time since the SN exploded. 
For f-b and b-f transitions, the rates are as in Eq.(18). For 
collisions, the rates are as in Eq.(19). 

4.3.2. Absorption of packets 

The absorption of r-packcts by lines is determined by the 
Sobolev optical depths given by Eq.(22). Absorption of 
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an r-packet to the continuum is determined by the con- 
ventional macroscopic absorption coefficient corrected for 
stimulated emission. 



4.3.3. Emission of packets 

The frequency of an emitted r-packet is decided as follows: 
for de-activation from a bound state i, the transition i — > j 
is selected with probability Aij (ijitiieij / E R , where E R is 
evaluated with Eq.(l) using the decay rates from Eq.(20), 
and the emitted packet is assigned frequency vT. - i.e., 
it is in the far red wing of a line whose emission profile 
is approximated by a delta function. The packet's next 
possible b-b transition is therefore with the next line to 
the redward of v^ (Abbott & Lucy 1985). 

For de-activation from a continuum state, the new fre- 
quency is, as in Sect. 4.2.3, selected by randomly sampling 
the energy distribution of the spontaneous recombination 
continua. 



4.3.4. Direction of propagation 

If an r-packet is emitted from a continuum state, the new 
direction of propagation is selected according to isotropic 
emission since the emission in this case is spontaneous. 
For de-activation from a bound state, the emission is 
also isotropic since, for homologous expansion, there is 
no kinematically-preferred direction. This is not true for 
a stellar wind. 



5. Convergence tests 

The Monte Carlo transition probabilities derived in Sect. 
3 are designed to reproduce asymptotically the emissiv- 
ity of an atomic species whose level populations are in 
statistical equilibrium. To test this, we now consider one- 
point problems with specified and fixed ambient condi- 
tions. Such tests sensibly precede application to a general 
NLTE problem, for then the local ambient conditions are 
everywhere being adjusted iteratively as the global solu- 
tion is sought. 

5.1. Fell 

In the initial tests, the Monte Carlo transition probabil- 
ities are applied to the model Fc II ion with N = 394 
levels used previously (Lucy 1999b) to investigate the ac- 
curacy of approximate treatments of line formation in 
SNe envelopes. The energy levels of the Fe II ion and 
the f-values for permitted transitions were extracted from 
the Kurucz-Bell (1995) compilation by M.Lennon (Mu- 
nich) . Einstein A- values for forbidden transitions are from 
Quinet et al.(1996) and Nussbaumer & Storey (1988). Col- 
lision strengths, needed for Sect. 5. 1.5, are from Zhang & 
Pradhan (1995) and van Regemorter (1962). 



5.1.1. Radiative excitation 

In the first Fe II test, we neglect collisional excitations 
and, as previously (Lucy 1999b), take the ambient ra- 
diation field determining the quantities Jj t in Eq.(20) 
to be WB„(T b ) with T b = 12500^ and dilution factor 
W = 0.5, corresponding to r = R. The density parameter 
is n(Fell) = 6.6 x 10 7 cto~ 3 , and the time since explosion 
is £e == 13 days. With parameters specified, this one-point 
statistical equilibrium problem - Eq.(4) for N — 1 levels 
plus a normalization constraint - is non-linear in the un- 
knowns rii because the rate coefficients in Eq.(20) depend 
on the rii through the Sobolev escape probabilities. Fortu- 
nately, repeated back substitutions give a highly accurate 

(x) 

solution n) in ~ 10 iterations. 



5.1.2. Monte Carlo experiment 



With nf' determined, the Fe II level emissivites E R and 
absorption rates A R can be computed from Eq.(l). We 
now test the Monte Carlo transition probabilities by see- 
ing how accurately they reproduce these values E R ; . Note 
that it is sufficient to test level emissivities since if these 
are exact so also are the line emissivities computed as de- 
scribed in Sect. 4.3.3. 

In the following Monte Carlo experiment, M packets of 
radiant energy are absorbed and subsequently emitted by 
a macro-atom representing a macroscopic volume element 
of Fe II ions in the ambient conditions specified above. 
The energies of these packets are taken to be equal and 
given by eo = A R /Af, where A R = J^. Af. The calculation 
proceeds step-by-step as follows: 

1) Mi = M A R j A R of the packets activate the macro- 
atom to internal state i. 

2) The transition probabilities p R , pi U and pu for a 
macro-atom in state i are computed from Eqs.(7) and (9). 

3) The transition probabilities sum to one, so each cor- 
responds to a segment (xk,Xk+i) of the interval (0,1). A 
particular transition is therefore selected by computing a 
random number x in (0,1) and finding in which segment 
it falls. 

4) If the selected transition is the de-activation of the 
macro-atom, we update E^ c to E^ +eo and then return 
to step 3) to process the next activation of state i, or to 
step 2) to process the first of the packets that activate the 
macro-atom to state i + 1. 

5) If the selected transition is an internal transition to 
state j, then we return to step 2) with j replacing i. 

6) When all M packets have been processed, the final 
elements of the vector Ef 10 are the estimates of the level 
emissivities E R 
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Fig. 2. Convergence test. The mean error 5 defined by 
Eq.(23) is plotted against TV, the number of packets in 
the Monte Carlo experiment. The open circles refer to the 
case where excitation energies are increased by 5eV. The 
straight line drawn by eye has slope = —0.5. Also indicated 
is the mean error when the level emissivities are assumed 
equal to the level absorption rates. 



5.1.3. Results of experiment 

As a single measure of the accuracy of the estimated level 
emissivities, we compute the quantity 



El*. 



MC 



E 



/£4' 



(23) 



This is the mean of the absolute fractional errors of the 
E^ c when weighted by Ef. 

Figure 2 shows the values of S, expressed as percentage 
errors, found in a series of trials with M increasing from 
10 4 to 10 7 . The values of 5 decrease monotonically with 
increasing A/", falling to 0.36 percent for M — 10 7 . More 
importantly, the errors accurately follow an A/" -1 ' 2 line, as 
expected if the only source of error are the sampling error 
at step 3) of the Monte Carlo experiments. Accordingly, to 
the accuracy of these experiments, macro-atoms obeying 
the transition probabilities derived in Sect. 3 do indeed 
reproduce the emissivity of a gas in statistical equilibrium. 

Also included in Fig. 2 are values of 6 obtained when 
the transition probabilities are computed with excitation 
energies ej increased by 5eV. This is to investigate the 
consequences of the dependence of the energy flow terms 
in Eq.(5) - and therefore also of the transition probabilities 
- on the zero point of the scale of excitation energy. These 
results also track an A/"" 1 / 2 line and so indicate that the 
predicted emissivities are asymptotically independent of 
the zero point. But since the open circles are marginally 



higher, there is an indication that increasing the zero point 
gives slighty less accurate emissivities at a given J\f. 

In the Monte Carlo codes for hot star winds and SNe 
cited in Sect. 1, line formation is treated approximately, 
with either resonant scattering or downward branching 
being assumed. For both assumptions, E^ = Af, corre- 
sponding to a macro-atom for which de-activation always 
immediately follows activation - i.e., pf = 1 for all i. In 
this case, as indicated on Fig. 2, 5 — 7.28 percent. Thus, 
when the points in Fig. 2 drop below this value, the suc- 
cess must be due to the internal, radiationless transitions 
governed by the probabilities pi U and pu. 

5.1.4. Distribution of jumps 

The above experiments show that despite the formidable 
complexity of its level structure the Fe II ion's repro- 
cessing of radiation is accurately simulated by the Monte 
Carlo transition probabilities. Nevertheless, from a com- 
putational standpoint, a remaining concern is how many 
internal transitions - or jumps - does this require? To an- 
swer this, the number of jumps before de-activation was 
recorded for each absorbed packet in the N — 10 7 trial 
and used to derive N(j), the number of packets requiring 
j jumps. 

From N(j), we find that the expected number of jumps 
is < j > = 2.19 and that the probability of immediate de- 
activation - i.e., zero jumps - is Pq = 0.425. Evidently, 
fears of numerous, time- consuming internal transitions 
are ill-founded. 

Figure 3 is a logarithmic plot of N(J). This reveals a 
power-law decline with increasing j but with alternating 
deviations indicating that an even number of jumps before 
de-activation is favoured. A simple model suggests the ori- 
gin of this curious behaviour. Consider a 3-level atom with 
£3 > £2 > ei = and suppose that level 2 is metastable 
with A21 = 0. Because B12 — 0, the macro-atom can only 
be activated to state 3; and because A21 = -B21 = 0, the 
macro-atom cannot de-activate from state 2. Moreover, 
since ei = 0, Eq.(9) gives P31 = P21 = 0, and so state 
1 of the macro-atom cannot be reached. Accordingly, fol- 
lowing activation at state 3, the macro-atom de-activates 
with probability p or jumps to state 2 with probability 

1 — p, from whence it returns to state 3 with probability 
P23 = 1. It is now simple to prove that the probabilty of j 
jumps before de-activation is Pj = p(l —p)^ 2 if j is even, 
and Pj — if j is odd. The Fe II ion's numerous low-lying 
metastable levels are presumably playing the role of level 

2 and thereby favouring an even number of jumps. 

Histograms N(j) have also been computed for two 
other cases. First, the above trial was repeated with the 
ej's increased by 5eV as in Sect. 5.1.3. This change in- 
creases < j > - to 4.54 - as expected since the probabilities 
Pi u and pn are thereby increased and pf" correspondingly 
decreased. Evidently, the standard choice of energy-level 
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Fig. 3. Histogram of N(j), the number of times in an ex- 
periment with J\f = 10 7 that the macro-atom underwent j 
internal transitions - or jumps - before de-activating with 
the emission of an energy packet. The mean number of 
jumps < j > and the probabilities of de-activation after 
j = — 2 jumps are indicated. 



zero point leads to the most computationally-efficient set 
of transition probabilities. 

In the second case, W is decreased from 0.5 to 0.067, 
corresponding to r = 2R. This change decreases < j > - 
from 2.19 to 1.29 - as expected given the weakening of the 
radiative excitation rates. 



5.1.5. Collisional excitation 

In the above experiment, the emission derives entirely 
from radiative excitation since collisions were neglected. 
Now we consider the opposite extreme by setting the am- 
bient radiation field to zero but including collisions. 

The only parameters of this test are the electron tem- 
perature and density, and these are assigned the values 
T e = 2 x 10 4 AT and N e = 10 8 cm" 3 . The resulting statis- 
tical equilibrium problem is linear and so solved without 
iteration. For this solution, accurate values of the level 
emissivities Ef- are again computed from Eq.(l). 

The next step is to derive estimates of the level 
emissivities by repeating the Monte Carlo experiment of 
Sect. 5. 1.2. The only changes needed are the following: 
first, since the solution has population inversions the gen- 
eral formulation of Sect. 4.1 must be adopted to avoid 
negative probabilities. 

Secondly, since a macro-atom is now always activated 
by a fc-packet, their energies are taken to be eo = A c /Af, 



where A c — J^-Af . Correspondingly, at step 1) of the 
experiment, A/i = JVAf'/A c . 

Thirdly, since a macro-atom can now de-activate by 
emitting cither an r- or a /c-packet, only the former re- 
sults in an updating of Ef lc . The emission of a /c-packet 
represents the return of energy eo to the therrmal pool. 

Apart from these changes, the convergence experiment 
proceeds as in Sects. 5.1.2 and 5.1.3. The result is a plot 
similar to Fig. 2, but with S = 0.19 percent for M = 10 7 . 
Evidently, the Monte Carlo transition probabilities are 
equally applicable to problems where collisional excitation 
is a source of emission. 

5.2. Hydrogen 

Although the Fe II experiments demonstrate the validity 
of the Monte Carlo transition probabilities, a test includ- 
ing b-f and f-b transitions is of interest. Accordingly, a con- 
vergence experiment at one point in a SN's envelope has 
also been carried out for a 15-lcvcl model of the H atom, 
with level 15 being the continuum k. The 14 bound levels 
correspond to principal quantum numbers n = 1 — 14, with 
each level having consolidated statistical weight g = 2n 2 . 
As for Fe II, the ambient radiation field incident on 
the blue wings of the b-b transitions is WB v (Tb), but 
now with Tf, = 6000AT and W — 0.067. However, beyond 
the Lyman limit, we assume zero intensity, so that pho- 
toionizations occur only from excited states. Correspond- 
ingly, recombinations to n = 1 are excluded on the as- 
sumption of immediate photoionization. Collisional exci- 
tations and ionizations are neglected. The density param- 
eter is N(H) = 1.88 x 10 9 cto~ 3 , the electron temperature 
T e = 4800AT, and the time since explosion £# = 10 days. 
With parameters specified, this non-linear statistical equi- 
librium problem can also be solved with repeated back 

(x) 

substitutions, giving a highly accurate solution n\ in 
~ 30 iterations. 

(x) 

With n\ determined, Monte Carlo experiments as de- 
scribed in Sect. 5.1.2 were carried out to test if level emis- 
sivities are also recovered in this case. In Fig. 4, two such 
trials, with Af = 10 4 and 10 5 , are compared with the ex- 
act solution. The results show that excellent agreement is 
achieved for AT = 10 5 . Note in particular the success with 
E^, which is the rate of ionization energy loss due to re- 
combinations, and with E^, whose very low value is due 
to the strong trapping of La photons. 



5.3. Alternative test of convergence 

Thus far, a Monte Carlo procedure has been used to vali- 
date the transition probabilities developed in Sect. 3. This 
has the advantage of following closely and therefore illus- 
trating their use in realistic NLTE calculations. But for 
feasible values of A/", sampling errors limit the accuracy of 
such tests. 
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6. Sensitivity tests 



The experiments of Sect. 5 demonstrate that, when COm- 
puted with the exact level populations n\ , the Monte 
Carlo transition probabilities applied to indivisible e- 
packets reproduce the exact level emissivities as Af — ► oo. 
But this success, though necessary, does not of itself im- 
ply that the technique will be successful when applied to 
NLTE problems. For example, if the Monte Carlo emissiv- 
ities were to undergo large changes in response to small 
changes in m, then we would reasonably suspect that the 
iterations inevitably required for a NLTE problem would 
converge very slowly - or might even diverge. On the other 
hand, if the emissivities are insensitive to changes in n,, 
then the prospects for successful applications are excel- 
lent. 



6.1. Fe II emissivities 



Fig. 4. Level emissivities (cgs) for Hydrogen. Results for 
trials with Af = 10 4 (large open circles) and Af = 10 5 
(small open circles) are compared with exact values (filled 
circles). The lines contributing to the level emissivities are 
indicated for n = 2 — 4. 

In order to test to higher precision, approximate level 
emissivities E\ m can be computed recursively according 
to the following scheme: 



4 (m) =pfE G 



(■'■■) 



(24) 



r=\ 



where pf is the radiative de-activation probability from 
Eq.(7) and G\ r is the increment at cycle r to the sum- 
mation approximating the rate at which level i gains en- 
ergy - i.e. the right-hand side of Eq.(5). This increment 
is derived from the previous increment by applying the 
transition probabilities from Eq.(9). Thus 



Gr ] = Y.p^ G . 



and the recursion cycles are initiated by setting 



(25) 



(26) 



This procedure is now applied to the Fe II test prob- 
lem of Sect. 5.1.1. As with that experiment, the accuracy 
of the vectors E\ m are measured by computing <5 defined 
by Eq.(23). For m = 17, 6 drops below the value 0.36 per- 
cent found in Sect. 5. 1.3 with Af = 10 7 - see Fig. 2. As the 
recursion procedure continues further, d decreases mono- 
tonically until at m ~ 60 it drops to a value of ~ 10~ 8 , 
at which point machine precision or accumulated roundoff 
errors halt further progress. This test clearly confirms and 
strengthens the earlier tests of the Monte Carlo transition 
probabilities. 



This crucial question of sensitivity can be investigated by 
repeating the calculations of Fe II emissivities reported in 

Sect. 5.1, but with m perturbed away from nf . A con- 
fa;) 
venicnt way of doing this is to replace n\ by the Boltz- 

mann distribution at excitation temperature T ex . Then, 

for given T ex , the corresponding level emissivities E^ 

are obtained from a Monte Carlo trial with Af = 5 x 10 6 

packets, and so are negligibly affected by sampling errors 

(cf. Fig.2) . 

Now, for the given T ex , we can also compute Ef", the 
level emissivities predicted by the fundamental formulae - 
Eqs.(l) and (20) in this case. This represents the standard 
approach to NLTE transfer problems whereby the radi- 
ation field is computed from the Radiative Transfer Eq. 
(RTE) with emissivity coefficients evaluated using the cur- 
rent estimates of m. Thus by comparing these two emis- 
sivity estimates E^ and E^ 1 , we can see whether this 
Monte Carlo technique is potentally capable of yielding a 
superior estimate of the radiation field. 

In Fig. 5, the quantities Ef 10 and Ef- obtained for 



r„ 



• (x) 

12500-fsT are plotted against E\ , the exact sta- 
tistical equilibrium level emissivities - i.e., the values cor- 
responding to nf . 

Remarkably, Figure 5 shows that the Monte Carlo 
emissivities are far less sensitive to the departure of m 

(x) 

from n\ than are the emissivities computed directly from 
the fundamental formula. For the most part, the Ef 10 are 
in error by < 0.1 dex, with little evidence of bias, while 
the Ef- are systematically offset by ~ +0.3 dex. 

To investigate whether this insensitivity is character- 
istic of the Monte Carlo procedure, the above test is now 
repeated with T ex ranging from 7500-ftT to 20000iiT and 
the resulting mean errors defined by Eq.(23) plotted in 
Fig. 6. We see that Ef gives reasonably accurate emissivi- 
ties only in the immediate neighbourhood of the minimum 
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Fig. 5. Sensitivity test. For a Boltzmann distribution over 
excited states at T ex = 12500i^T, the level cmissivities (cgs) 
obtained with the Monte Carlo transition probabilities 
(filled circles) and with the basic formula (open circles) 
are plotted against the exact emissivities obtained with 

(x) 

n\ . The Monte Carlo emissivities derive from a trial with 
M = 5 x 10 6 packets. 



at T ex ~ 11250AT. On the other hand, the values Ef c are 
accurate to SS 0.1 dex across the entire range. 

The causes of these astonishing differences in sensitiv- 
ity are of considerable interest. For Ef, the strong sen- 
sitivity to T ex is readily understood. Because the sum 
Ru^ie oc rii, an error in the population of the emitting 
level translates directly into an error in Ef- . 

Now consider E^ c . This quantity is determined by the 
rate at which active macro- atoms reach state i, and this 
happens by direct absorptions of packets into this state or 
by transitions from other states. Either way, the accuracy 
of the source vectors Af and A^ is clearly fundamental 
to the accuracy of the vector Ef 10 . But the dominant 
contributors to the elements of these source vectors - see 
Eqs.(l) and (2) - are transitions from the ground state 
and from low-lying metastable levels, and the estimated 
populations of these levels are unlikely to be seriously in 
error. In particular, with an assumed Boltzmann distri- 
bution over excited states, the m of these low levels is 



insensitive to T PT and do not differ much from 



» 



In 



contrast, the populations of high levels arc quite likely to 
be badly estimated and are acutely sensitive to T ex . 

6.2. Comments 

Another way of appreciating the differences in these ap- 
proaches to calculating cmissivities is as follows. The 
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Fig. 6. Sensitivity test. Logarithmic errors of the emissiv- 
ity vectors Ef IC and Ef evaluated for Boltzmann distri- 
butions over excited states plotted against T ex . The Monte 
Carlo emissivities derive from trials with N = 5 x 10 6 
packets. 



Monte Carlo procedure applies only to a state of statistical 
equilibrium and, as such, constrains every level's emissiv- 
ity to be consistent with the rates of processes populating 
that level. In contrast, the fundamental emissivity formula 
applies also to states out of statistical equilibrium and so 
takes no account of whether the levels' populations can 
be maintained. Accordingly, with this Monte Carlo tech- 
nique, the principle of statistical equilibrium is incorpo- 
rated (approximately) as the radiation field is being calcu- 
lated. On the other hand, when emissivities are computed 
from the fundamental formula, any consideration of sta- 
tistical equilibrium is effectively being deferred until the 
updated radiation field has been determined. 

The likely beneficial impact of this insensitivity on the 
iterations needed to derive NLTE solutions is worth stress- 
ing. With the conventional RTE approach, an erroneously 
overpopulated upper level i pollutes the radiation field 
with spurious line photons at frequencies i/jj (j < i), and 
these are sources of excitation for level i when level popu- 
lations are next solved for. Similarly, an erroneously over- 
populated upper ion pollutes the radiation field with re- 
combination photons that are subsequent sources of pho- 
toionization for the lower ion. To some degree, therefore, 
such errors are self-perpetuating and so are not rapidly 
eliminated. This persistency contributes to the slow con- 
vergence typical of NLTE codes. In contrast, with the 
Monte Carlo approach, this pollution does not happen and 
so - for sufficiently large N - a high quality radiation field 
is obtained immediately provided that the initial popula- 
tions of the low-lying levels are estimated sensibly 
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7. Implementation 

The Monte Carlo transition probabilities allow statistical 
equilibrium to be incorporated into the calculation of radi- 
ation fields for NLTE problems. Moreover, this is achieved 
without imposing the constraint of radiative equilibrium. 
Accordingly, in principle at least, the technique applies 
equally to problems with non-radiative heating, such as 
stellar chromospheres. 

7.1. Radiative equilibrium 

In the absence of non-radiative heating, a NLTE transfer 
problem must be solved subject to the constraint of ra- 
diative equilibrium. The incorporation of this additional 
constraint into the macro-atom formalism is readily un- 
derstood. First suppose that collisional processes are ne- 
glected. The absorbed and the emitted e-packets are then 
always r-packcts and they have identical energies - see 
Fig.l. Thus, the constraint of radiative equilibrium is 
obeyed rigorously since it holds exactly for every activa- 
tion - de-activation event, all of which are of the form 
r — > A* — > r, where A* denotes an active macro-atom. 
Note also that since active macro-atoms do not appear 
spontaneously within the computational domain (D), ev- 
ery Monte Carlo quantum's interaction history starts and 
ends as an r-packet crossing a boundary of D. 

Now suppose that collisions are included. In this case, 
a macro-atom activated by an r-packet can de-activate 
itself by emitting a fc-packct, so that radiative equilib- 
rium no longer holds exactly for each individual activa- 
tion - de-activation event. However, the emitted fc-packet 
is re-absorbed in situ by another macro-atom and thereby 
(eventually) converted into an r-packet. Since this has the 
same energy as the original r-packet, radiative equilib- 
rium holds for every sequence of in situ events that starts 
with the absorption of an r-packet and ends with the 
next emission of an r-packet. A typical in situ sequence is 
r — > A* — > k — > A* — > k — > A* — ► r. If such sequences are 
abbreviated as r — > [A*] — > r, we see that the inclusion 
of collisions has not fundamentally changed the procedure 
and that radiative equilibrium is still rigorously obeyed. 

7.2. Non-radiative heating 

In the presence of non-radiative heating, the NLTE prob- 
lem is not subject to the additional constraint of radiative 
equilibrium. Statistical equilibrium is incorporated with 
the macro-atom formalism as before, and the challenge 
now is to incorporate the creation of radiant energy within 
D due to the additional heating. This is accomplished 
by allowing for the spontaneous and random appearance 
within D of active macro-atoms with their number, loca- 
tions and internal states i all controlled by the collision 
source vector Af - cf. Sect. 5. 1.5. Note that because this 
sampling of Af takes full account of the collisional cre- 



ation of excitation, the emission of a fc-packet is not now 
followed by its in situ re-absorption; instead, the interac- 
tion history of that Monte Carlo quantum then ends and 
its energy is added to the thermal pool (cf. Sect. 5. 1.5.). 
The radiation field generated by this procedure is not 
divergence-free but reflects the collisional creation of ra- 
diant energy due to an elevated temperature profile main- 
tained by the non-radiative heating. 

8. Conclusion 

The limited aim of this paper has been to see if Monte 
Carlo transfer codes whose quanta are indestructable en- 
ergy packets can be constructed without resorting to sim- 
plified treatments of line formation. To this end, the con- 
cept of a macro-atom has been introduced and rules estab- 
lished governing its activation and de-activation as well as 
its transitions between internal states. These rules - the 
Monte Carlo transition probabilities - have been derived 
by demanding that the macro-atom's emission of r-packcts 
asymptotically reproduces the local emissivity of a gas in 
statistical equilibrium; and these rules' validity has been 
confirmed with one-point test problems. 

Evidently, the next step is to implement these tran- 
sition probabilities in a code to solve a realistic NLTE 
problem for a stratified medium and thus to investigate 
the practicality of this technique for problems of current 
interest. In a companion paper, a Monte Carlo NLTE code 
treating the formation of H lines in a Type II SN envelope 
will be described and used to illustrate the convergence 
behaviour of iterations to obtain both the level popula- 
tions and the temperature stratification. 
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